module myprecision 
 implicit none  
    character(len=20)          :: lt_my_version = "V45        "

    !integer, parameter         :: sp = KIND(1.0)
    integer, parameter         :: dl = KIND(1.d0)
    integer, parameter         :: pr = KIND(1.d0)
    real(pr), parameter        :: fourPI=4.0_pr*3.1415926535897932384626433_pr    
    real(pr), parameter        :: myee = 2.71828183_pr
    real(pr), parameter        :: C_light = 2.99792458e8_pr
    real(pr), parameter        :: Cfrom1e = 1.60217733e-19_pr  !库仑电荷
    real(pr), parameter        :: Jfrom1eV = 1.60217653e-19_pr ! 1eV = ? J(焦耳)
    real(pr), parameter        :: eVfrom1Ry = 13.606_pr 
    real(pr), parameter        :: epsilon0 = 8.854187817e-12_pr
    real(pr), parameter        :: one_au = 0.5291772108e-10_pr    ! 1a.u.=? m ,bohr半径
    real(pr), parameter        :: re0 = 2.8179402894e-15_pr     ! 电子经典半径
    real(pr), parameter        :: alpha_jingxi = 137.035999074_pr !精细结构常数
end module myprecision 



module module_mychi2 
 use  myprecision
 implicit none  

    !character(len=100)          :: fname_exp = "si.dat"  
    character(len=100)          :: fname_exp = "si.Spe"  
    integer                     :: nlines_exp_ini
    real(pr)                    :: dt_exp  !! ns  道宽
    real(pr)                    :: line3_exp
    real(pr)                    :: fbl_exp      !! ns 分辨录
    integer,allocatable         :: ltpu_exp_ini(:) 


    integer,parameter           :: nleft =  5
    integer,parameter           :: nright = 800 
    integer,parameter           :: nleft_theoini = 100
    real(pr)                    :: ltpu_exp(1:nleft+nright+1) 
    real(pr)                    :: ltpu_theo(1:nleft+nright+1) 
    real(pr)                    :: ltpu_d_exp_theo(1:nleft+nright+1) 
    integer,parameter           :: nbin_perexp= 6 !! 实验数据每道被分割为更多,为了更精细的理论谱  
    real(pr)                    :: ltpu_background
    integer                     :: dn_background = 150 

    real(pr),allocatable        :: gauss_array(:)
    integer                     :: gauss_nleft
    integer                     :: gauss_nright
 
    integer,parameter           :: ntau = 3
    real(pr)                    :: value_tau(1:ntau) != (/6._pr, 0.385_pr, 0.220_pr /)
    real(pr)                    :: intensity(1:ntau) != (/0.005_pr , 0.10_pr, 0._pr /)


    TYPE                        :: type_interplate
        integer                 :: nmin, nmax
        real(pr)                :: xmin, xmax, xdelta
        real(pr)                :: resolution  !! unit : ps
        real(pr),allocatable    :: array(:)
        real(pr),allocatable    :: array_c(:)  
        real(pr),allocatable    :: array_fitexp(:)  
        !! : after convolution of array and gauss function
    end TYPE 
    TYPE(type_interplate)    :: dNofdt !! 光子计数（T） !! /ps
    


 contains 


    subroutine countlines_v2(fname,nlines)
       character(len=100),intent(in)     :: fname
       integer,intent(out)                 :: nlines
       character(len=100)                :: temp
       nlines = 0
       open(1,file=fname) 
       do while(1==1) 
           read(1,*,err=11,end=12) temp
           nlines = nlines + 1
   !        write(*,*) temp
       enddo
11     write(*,*) " error occurs when calling countlines";    stop
12     close(1);return !;write(*,*) "nlines=", nlines 
    end subroutine


    subroutine read_expdata
        character(len=10)        :: temp_10c
        integer                 :: i 

        allocate(ltpu_exp_ini(1:nlines_exp_ini-4 ) )
        open(111,file=trim(adjustl( fname_exp ) ) )
            read ( 111, * )  temp_10c
            read ( 111, * )  temp_10c  !dt_exp
            read ( 111, * )  temp_10c  !line3_exp
            read ( 111, * )  temp_10c  !fbl_exp
            do i = 1, size(ltpu_exp_ini)
                read ( 111, FMT=*,err=111,end=112)  ltpu_exp_ini(i)
            enddo 
        close(111)
        !print*, fbl_exp 
        return 

111  print*, "   error occurs when 11111461663 " ;stop
112  print*, " end 49813613";close(123) ;return !;write(*,*) "nlines=", nlines 
    end subroutine 


    subroutine mod_data
        integer             :: i 
        integer             :: i_maxloc_exp
        character(len=100)  :: temp_100c
        !
        call countlines_v2(fname_exp,nlines_exp_ini)
        print*, "nlines_exp_ini = ",  nlines_exp_ini 
        call read_expdata
        i_maxloc_exp = maxloc(ltpu_exp_ini,1 )
        !
        do i = 1, nleft+nright+1
            ltpu_exp(i) = ltpu_exp_ini(i_maxloc_exp-nleft-1+i) 
        enddo

        ltpu_background = sum(ltpu_exp_ini( i_maxloc_exp+nright-dn_background : & 
                                        i_maxloc_exp+nright+dn_background) ) /dn_background/2._pr
        print*, "ltpu_background = ", ltpu_background

    end subroutine 
    


    subroutine get_ltpu_theo_ini(temp_ltpu_theo)
        real(pr),intent(inout)  :: temp_ltpu_theo(1: (nleft+nright+1+nleft_theoini)*nbin_perexp )  
        integer                 :: i , i_ntau
        do i_ntau = 1, ntau
        do i = (nleft_theoini*nbin_perexp+1), size(temp_ltpu_theo)
            temp_ltpu_theo(i) = temp_ltpu_theo(i) +  intensity(i_ntau)/value_tau(i_ntau) &
                             / exp( ( i-(nleft_theoini*nbin_perexp+1) )  &
                             * dNofdt%xdelta / value_tau(i_ntau) )
        enddo            
        enddo

    end subroutine



    function GaussSimpson(AS,BS,NS,FWHM) 
       ! Integral of Gauss_Func form AS to BS  
       IMPLICIT NONE
       real(pr), EXTERNAL             :: FX
       real(pr)                     :: FWHM 
       real(pr)                     :: GaussSimpson,AS,BS
       real(pr)                     :: HS
       real(pr)                     :: zS1,zS2
       real(pr)                     :: SS,fS1,fS2
       INTEGER                     :: NS,I
   
       HS=(BS-AS)/NS
       zS1=AS
       zS2=zS1+HS
       SS=0.0_pr
       fS1=Gauss_Func(zS1,FWHM)
       fS2=Gauss_Func(zS2,FWHM)
       !N=40
       SS = HS*(fS1+4.0_pr*(Gauss_Func( (zS1+zS2)*0.5_pr,FWHM) ) &    
                               + fS2)/6.0_pr
       DO I = 2, NS
           zS1=zS2
           fS1=fS2
           zS2=zS2+HS
           fS2=Gauss_Func(zS2,FWHM)
           SS = SS + HS*(fS1+4.0_pr*Gauss_Func( (zS1+zS2)*0.5_pr,FWHM) &
                           +fS2)/6.0_pr
       END DO
       GaussSimpson = SS
       RETURN
    end function
   
   
    function Fx_3points(x,x1,F1,x2,F2,x3,F3)
       real(pr)                :: x,x1,F1,x2,F2,x3,F3
       real(pr)                :: Fx_3points
       Fx_3points = (F1*(x - x2)*(x - x3))/((x1 - x2)*(x1 - x3)) & 
                   +(F2*(x - x1)*(x - x3))/((x2 - x1)*(x2 - x3)) &
                   +(F3*(x - x1)*(x - x2))/((x3 - x1)*(x3 - x2))
    end function


    !!!!!!!!!!!!!!!!!!!!!!!!!!原始的dN/dt 的插值（没有卷积）!!!!!!!!!!
    function Intplt_array(x,vars)
       IMPLICIT NONE
       TYPE(type_interplate)                 :: vars 
       real(pr)                            :: x
       real(pr)                            :: Intplt_array
       integer                                :: n_x
       real(pr)                            :: x_n,x_pre, x_next
       !allocate(vars%array(vars%nmin:vars%nmax) )
       !deltax = (vars%xmax-vars%xmin)/(vars%nmax-vars%nmin)
       if ( x .LT. vars%xmin .or. x .GT. vars%xmax) then 
           write(*,*) "error No.1234056000999 ";stop
           return 
       elseif ( x .GT. vars%xmin+vars%xdelta  & 
               .and. x .LT. vars%xmax-vars%xdelta ) then
           n_x = int( (x-vars%xmin)/vars%xdelta + 0.5_pr) !+ vars%nmin
           x_n = n_x * vars%xdelta + vars%xmin
           x_pre = x_n - vars%xdelta 
           x_next = x_n + vars%xdelta
       elseif (x .LE. vars%xmin+vars%xdelta) then
           n_x = vars%nmin + 1 
           x_n = vars%xmin + vars%xdelta !* n_x   
           x_pre = x_n - vars%xdelta 
           x_next = x_n + vars%xdelta
       else!if (x .LT. xmax-deltax)
           n_x = vars%nmax - 1 
           x_n = vars%xmax - vars%xdelta !<-n_x * vars%xdelta + vars%xmin 
           x_pre = x_n - vars%xdelta 
           x_next = x_n + vars%xdelta
       endif
       Intplt_array = Fx_3points(x,x_pre,vars%array(n_x-1), &
                   x_n,vars%array(n_x),x_next,vars%array(n_x+1))
    end function 
    
    
    function gauss_func(x,fwhm) !!! one dimention 
       IMPLICIT NONE
       real(pr)        :: Gauss_Func,x
       real(pr)        :: FWHM !! 半高全宽 
       real(pr)        :: sigma
       sigma  = FWHM/2.354820047_pr != FWHM/(2._pr*1.1774100236_pr)
       !Gauss_Func = exp(-0.5_pr*x**2.0_pr/sigma**2.0_pr)/sqrt(fourPI/2.0_pr*sigma) 
       Gauss_Func = exp(-0.5_pr*x**2.0_pr/sigma**2.0_pr) &
                       / (sqrt(fourPI/2.0_pr)*sigma )
    end function
    
    
    function Convolution(Ntimes,Gauss,FWHM,vars_intplt,x,vars)
       IMPLICIT NONE
       real(pr)                    :: Convolution
       integer                        :: Ntimes
       real(pr),external            :: Gauss
       real(pr)                    :: FWHM 
       real(pr),external            :: vars_intplt 
       TYPE(type_interplate)        :: vars
       real(pr)                     :: x!, alpha
       real(pr)                     :: AS, BS, HS, normalization
       real(pr)                     :: a1,a2,a_half ,xa1,xa2,xa_half
       real(pr)                     :: f1,f2,f_half
       INTEGER                     :: I
       !Fx_star(alpha) = vars_intplt(x - alpha,vars) * gauss_func(alpha,FWHM) 
       AS = vars%xmin
       BS = vars%xmax
    
       if (1==0) then
           if (int( (x-vars%xmin)/vars%xdelta + 0.5_pr) &
               == (int( (vars%nmax-vars%nmin)/2._pr) + vars%nmin) ) then
               !== 0 ) then
               normalization = GaussSimpson(AS,BS,Ntimes,FWHM) 
               if (abs(normalization-1._pr) .GE. 1.E-9_pr) then
               print*," Err NO.290889, Gauss Normalization :",normalization
               endif
           endif
       endif
       !normalization = 1.0_pr 
    
       Convolution = 0.0_pr 
       HS = (BS-AS)/Ntimes 
       a1 = vars%xmin 
       a2 = a1 + HS 
       a_half = (a1 + a2)*0.5_pr 
       xa1 = x - a1 
       xa2 = x - a2
       xa_half = x - a_half 
       f1 = vars_intplt(a1,vars) * Gauss(xa1,FWHM) !/ normalization
       f2 = vars_intplt(a2,vars) * Gauss(xa2,FWHM) !/ normalization
       f_half = vars_intplt(a_half,vars) & 
                   * Gauss(xa_half,FWHM)        !/ normalization
       Convolution = Convolution + HS*(f1+4.0_pr*f_half+f2)/6.0_pr
    
    
       DO I = 2, Ntimes
           a1 = a2
           a2 = a1 + HS 
           a_half = (a1 + a2)*0.5_pr 
           xa1 = x - a1 
           xa2 = x - a2 
           xa_half = x - a_half 
    
           f1=f2
           f2 = vars_intplt(a2,vars) * Gauss(xa2,FWHM) !/ normalization
           f_half = vars_intplt(a_half,vars) &
                       * Gauss(xa_half,FWHM)           !/ normalization 
           Convolution = Convolution + HS*(f1+4.0_pr*f_half+f2)/6.0_pr
       END DO
       RETURN
    end function 



    function get_gauss_array(myFWHM,myxdelta,mynleft,mynright)
        real(pr),intent(in)       :: myFWHM
        real(pr),intent(in)       :: myxdelta
        integer,intent(in)        :: mynleft, mynright
        real(pr),allocatable      :: get_gauss_array(:) 
        real(pr)                  :: myAS, myBS
        integer                   :: i 

        allocate(get_gauss_array(-mynleft:mynright) )
        do i = -mynleft , mynright   
            myAS = i*myxdelta - myxdelta/2._pr 
            myBS = i*myxdelta + myxdelta/2._pr  
            get_gauss_array(i) = GaussSimpson(myAS,myBS,8,myFWHM) 
        enddo
    
    end function 

    function mychi2(p1_dt_exp,p2_fbl_exp ,p3_tau1,p4_intsty1  ,p5_tau2,p6_intsty2  ,p7_tau3 )
        real                        :: p1_dt_exp
        real                        :: p2_fbl_exp
        real                        :: p3_tau1
        real                        :: p4_intsty1
        real                        :: p5_tau2
        real                        :: p6_intsty2
        real                        :: p7_tau3
        real                        :: p8_intsty3
        real                        :: p9_tau4
        real                        :: p10_intsty4
        !
        real(pr)                    :: mychi2
        real(pr)                    :: mychi2_best
        integer                     :: i_ntau
        integer                     :: i
        integer                     :: id  !! for dNofdt , WofLT
        real(pr)                    :: T_id
        integer                     :: imaxloc_fitexp  !理论谱粗化后的最大数指标 dNofdt%array_fitexp
        integer                     :: imaxloc_exp  !实验谱最大数指标 ltpu_exp
        integer                     :: sum_imaxloc_fitexp  !理论谱粗化后的最大数指标 dNofdt%array_fitexp
        integer                     :: sum_imaxloc_exp  !实验谱最大数指标 ltpu_exp
        integer                     :: nsum_imaxloc     !理论实验谱的最大位置数据规约时所选取的数据量
        integer                     :: nmax_chi2   !  计算chi2的最高道数
        integer                     :: nmin_chi2   !  计算chi2的最低道数
        integer                     :: nshift_theo_ini  ! 粗细理论谱之间变换时的偏移量


        value_tau(1:ntau)  = (/p3_tau1 ,p5_tau2 ,p7_tau3 /) !p_tau(1:ntau)   !(/6._pr, 0.385_pr, 0.220_pr /)
        intensity(1:ntau-1)= (/p4_intsty1 ,p6_intsty2 /)   !p_intsty(1:ntau-1) !(/0.005_pr , 0.10_pr, 0._pr /)
        !value_tau(1:ntau) = (/6._pr,     0.385_pr , 0.220_pr /)
        !intensity(1:ntau) = (/0.01_pr , 0.10_pr  ,  0.89_pr /)
        !
        intensity(ntau) = 1._pr 
        do i_ntau = 1, ntau-1
            intensity(ntau) = intensity(ntau) - intensity(i_ntau) 
        enddo
        if ( abs(sum(intensity) - 1._pr) > 1.E-7_pr )   stop 
        !write(*,'(A,3G13.5)'), " value_tau = ", value_tau
        !write(*,'(A,3G13.5)'), " intensity = ", intensity


        !print*, "nleft, nright ,nleft_theoini = ", nleft, nright, nleft_theoini 
        !print*, "nbin_perexp   = ", nbin_perexp
        if( .not. allocated(ltpu_exp_ini) ) call mod_data 
        !open(112,file="out_ltpu_exp_mod.txt")
        !    write(112,'(1E15.7)') (ltpu_exp(i),i= 1, (nleft+nright+1)  ) 
        !close(112)
        !open(112,file="out_ltpu_exp_mod_-bground.txt")
        !    write(112,'(1E15.7)') (ltpu_exp(i) - ltpu_background ,i= 1, (nleft+nright+1)  ) 
        !close(112)
        dt_exp  = p1_dt_exp 
        fbl_exp = p2_fbl_exp
        !print*, "dt_exp    = ", dt_exp    
        !print*, "fbl_exp   = ", fbl_exp    



        dNofdt%nmin = 1   !!!  因该不能改变 
        dNofdt%nmax = (nleft+nright+1+nleft_theoini)*nbin_perexp  !!在理论谱零点左边另取nleft_theoini个道
        dNofdt%xdelta = dt_exp/nbin_perexp 
        dNofdt%xmin = 0.0_pr   !! ns 
        dNofdt%xmax = (dNofdt%nmax-dNofdt%nmin)*dNofdt%xdelta + dNofdt%xmin    !! ns 
        dNofdt%resolution = fbl_exp   
        !
        allocate(dNofdt%array(dNofdt%nmin:dNofdt%nmax) )
        dNofdt%array = 0.0_pr
        allocate(dNofdt%array_c(dNofdt%nmin:dNofdt%nmax) )
        dNofdt%array_c = 0.0_pr
        allocate(dNofdt%array_fitexp( 1:  (nleft+nright+1+nleft_theoini) ) )
        call get_ltpu_theo_ini(dNofdt%array)
        !open(112,file="out_ltpu_theo_ini.txt" )
        !    write(112,'(1E15.7)') (dNofdt%array(i),i=1, (nleft+nright+1+nleft_theoini)*nbin_perexp ) 
        !close(112)


        gauss_nleft = int(5*dNofdt%resolution/dNofdt%xdelta)
        gauss_nright = gauss_nleft 
        allocate( gauss_array( -gauss_nleft : gauss_nright ) )
        gauss_array = get_gauss_array(dNofdt%resolution,dNofdt%xdelta, gauss_nleft, gauss_nright)
        !open(112,file="out_gauss.txt")
        !    write(112,'(1E15.7)') (gauss_array(i),i= -gauss_nleft, gauss_nright ) 
        !close(112)
        !
        dNofdt%array_c = 0._pr
        do id = dNofdt%nmin, dNofdt%nmax
            do i = id-gauss_nleft, id+gauss_nright
                if (  dNofdt%nmin <= i <= dNofdt%nmax ) then 
                    dNofdt%array_c(id) = dNofdt%array_c(id) + dNofdt%array(i) * gauss_array(i-id)
                endif
            enddo
        enddo 
        !open(112,file="out_ltpu_theo_conv.txt")
        !    write(112,'(1E15.7)') (dNofdt%array_c(i),i=1, (nleft+nright+1+nleft_theoini)*nbin_perexp ) 
        !close(112)


        mychi2_best = -1.0_pr
      do nshift_theo_ini = 0, nbin_perexp-1
        !print*, cshift( (/1,2,3,4/), -1)
        dNofdt%array_c = cshift(dNofdt%array_c, -1 ) 
        dNofdt%array_fitexp = 0.0_pr 
        do id = 1, (nleft+nright+1+nleft_theoini) 
            do i = (nbin_perexp*(id-1)+1), nbin_perexp*(id) 
                 dNofdt%array_fitexp(id) = dNofdt%array_fitexp(id) + dNofdt%array_c(i) 
            enddo
        enddo
        
         
        imaxloc_fitexp      = maxloc(dNofdt%array_fitexp,1)
        imaxloc_exp         = maxloc(ltpu_exp,1)
        ! 利用最大值附近的数据来规约 理论谱与实验谱
        !nsum_imaxloc        = 5
        !sum_imaxloc_fitexp  = sum(dNofdt%array_fitexp(imaxloc_fitexp-nsum_imaxloc: & 
        !                        imaxloc_fitexp+nsum_imaxloc) ) 
        !sum_imaxloc_exp     = sum(ltpu_exp(imaxloc_exp-nsum_imaxloc:imaxloc_exp+nsum_imaxloc) ) 
        !dNofdt%array_fitexp = dNofdt%array_fitexp / sum_imaxloc_fitexp * sum_imaxloc_exp  
        do i = 1, (nleft+nright+1) 
            ltpu_theo(i) = dNofdt%array_fitexp(imaxloc_fitexp-imaxloc_exp+i)  
        enddo 
        ltpu_theo = ltpu_theo/sum(ltpu_theo)* sum(ltpu_exp-ltpu_background)
        ltpu_d_exp_theo = ltpu_theo - ltpu_exp + ltpu_background 


        !nmax_chi2 = nleft + nright +1
        !nmax_chi2=100
        !do while ( ltpu_theo(nmax_chi2) > ltpu_background )
        !    nmax_chi2 = nmax_chi2 + 1 
        !enddo 
        nmax_chi2 = 501 !!int(20._pr/dt_exp)  ! 约 XXX ns之后不计
        ! 
        nmin_chi2 = 2   
        !print*, "nmin_chi2 = ", nmin_chi2 
        !print*, "nmax_chi2 = ", nmax_chi2 
        !print*, "size(ltpu_exp(nmin_chi2:nmax_chi2) ) ", size(ltpu_exp(nmin_chi2:nmax_chi2)) 
        mychi2 = sum( ltpu_d_exp_theo(nmin_chi2:nmax_chi2)**2 / ltpu_exp(nmin_chi2:nmax_chi2) )
        !!mychi2 = sum( (ltpu_exp(nmin_chi2:nmax_chi2)-ltpu_background) **2 / ltpu_background  )
        
        
        if ( mychi2_best > mychi2 .or. mychi2_best == -1.0_pr  ) then 
            mychi2_best = mychi2 
            !open(112,file="out_ltpu_theo_final.txt")
            !    write(112,'(1E15.7)') (ltpu_theo(i),i= 1, (nleft+nright+1)  ) 
            !close(112)
        endif  
      enddo
        mychi2=mychi2_best
        write(*,'( 2E15.7 , 6E13.5 )') mychi2 / size(ltpu_exp(nmin_chi2:nmax_chi2)), p1_dt_exp,p2_fbl_exp &
                        ,p3_tau1 ,p4_intsty1   ,p5_tau2 ,p6_intsty2   ,p7_tau3 
        !print*,  "mychi2_best = ", mychi2_best 
        !print*,  "mychi2/dof  = ", mychi2_best / size(ltpu_exp(nmin_chi2:nmax_chi2)) !/ (nleft+nright+1)


        deallocate(dNofdt%array )
        deallocate(dNofdt%array_c )
        deallocate(dNofdt%array_fitexp )
        if( allocated(gauss_array) )  deallocate(gauss_array)
    end function

end module module_mychi2



program main 
 use  myprecision
 use  module_mychi2
 implicit none
 real(pr)            :: main_tmp, para1

 main_tmp = mychi2(0.03  ,0.28   ,3. ,0.005   ,0.385 ,0.15  ,0.220 )
 main_tmp = mychi2(0.025 ,0.25   ,2. ,0.005   ,0.385 ,0.17  ,0.220 )

end program 

